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We study the thermodynamics of clean, layered superconductor/ferromagnet nanostructures using 
fully self consistent methods to solve the microscopic Bogoliubov-deGennes equations. From these 
self-consistent solutions the condensation free energies are obtained. The trilayer SFS junction 
is studied in particular detail: first order transitions between and 7r states as a function of the 
temperature T are located by finding where the free energies of the two phases cross. The occurrence 
of these transitions is mapped as a function of the thickness dp of the F layer and of the Fermi 
wavevector mismatch parameter A. Similar first order transitions are found for systems with a larger 
number of layers: examples are given in the 7 layer (3 junction) case. The latent heats associated 
with these phase transitions are evaluated and found to be experimentally accessible. The transition 
temperature to the normal state is calculated from the linearized Bogoliubov-deGennes equations 
and found to be in good agreement with experiment. Thus, the whole three dimensional phase 
diagram in T, dp, A space can be found. The first order transitions are associated with dips in the 
transition temperature T c to the non-superconducting state, which should facilitate locating them. 
Results are given also for the magnetic moment and the local density of states (DOS) at the first 
order transition. 

PACS numbers: 74.45.+C, 74.25.Bt, 74.78.Fk 



I. INTRODUCTION 

The investigation of systems involving ferromagnet (F) 
and superconductor (S) junctions is an active compo- 
nent of superconductor-based spintronics 1 research. A 
broad array of interesting effects arises in S /F nanostruc- 
tures, which opens doors for nanotechnologies and associ- 
ated devices and applications that may offer benefits be- 
yond current superconducting devices such as standard 
Josephson junctions. Advances in fabrication techniques 
permit growth of ferromagnet and superconductor layers 
in the form of junctions and heterostructures smooth up 
to the atomic scale. 

The arrangement of consecutive F and S layers, as in 
SFS junctions, results in competition between magnetic 
and superconducting orderings. Superconducting corre- 
lations can leak into the ferromagnet while spin polariza- 
tion can extend into the superconductor: these are the 
now well established S/F proximity effects. 2 ' 3 The phase 
coherence embodied in the superconducting correlations 
becomes modified in the F regions. The exchange energy 
in the ferromagnet shifts the kinetic energies of the quasi- 
particles constituting the Cooper pairs and subsequently 
a new superconducting state arises whereby the center of 
mass momentum of the pair is nonzero. 4 This results in 
a spatially decaying pair amplitude that oscillates over 
a characteristic length scale much smaller than the su- 
perconducting coherence length. The modulating pair 
amplitude within the magnet indirectly links adjacent S 
layers, and thus proximity effects in F cause local os- 
cillations in physically relevant single-particle quantities, 
including the magnetization 5,6 and density of states 7,8 



(DOS). Similarly, in the S material the magnet locally 
polarizes the superconductor, causing a monotonic de- 
cline in the pairing correlations near the interface over 
an extended region. The associated spin-split Andreev 
quasiparticle states also lead to interesting local behavior 
in the DOS and magnetic moment in the superconduc- 
tor. The nontrivial behavior of the proximity effects in 
these structures plays a central role in the competition 
between the magnetic and superconducting order. 

The modification of the superconducting phase coher- 
ence due to proximity effects in clean multilayers consist- 
ing of one or more successive SFS junctions is particu- 
larly striking. On the atomic level, the pair amplitude 
is a smoothly varying function of the spatial coordinates. 
Depending on the values of certain parameters (such as F 
layer width, dp) the damped oscillatory pair amplitude 
in the F layer may arrange itself in such a manner that is 
energetically favorable for its sign to change from one of 
the S layers to the next, yielding a so-called 7r-junction, 
as first proposed long ago. 9 If the pair amplitude does 
not change sign between S layers, it is an ordinary or 
0-junction. There is a rich and broad parameter space 
that then enables a certain level of control over the com- 
peting magnetic and superconducting orderings, allow- 
ing one to increase or diminish the proximity effects that 
dominate the relative SFS coupling. The actual equilib- 
rium state (0 or 7r) is dependent upon several variables, 
including predominantly the F region's material charac- 
teristics and the temperature, T, all of which ultimately 
determine the pair amplitude modulation in the magnet. 
A system comprised of a larger number of SFS sequences 
results in a greater number of possible or n junction 
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combinations. 

The transitions between and 7r states can be explored 
through the signatures of a variety of physical parame- 
ters. Experimental study of this question has focused 
primarily on measurements of the critical current I c 10-18 
and, thermodynamically, on the critical temperature 19-24 
of the transition to the normal state, T c . Evidence of 
«-> 7r transitions can be seen in the SFS Josephson 
coupling, which manifests itself in the vanishing of I c , 
although higher order harmonics in the current-phase 
relationship can modify this. 25 Measurements 19-24 as a 
function of d,F have shown that T c , which is of course 
smaller than T°, the critical temperature for bulk S ma- 
terial, oscillates as a function of (If, confirming theoreti- 
cal predictions 26,27 based upon the semi-classical Usadel 
equations. Intrinsically linked to this phenomenon are 
damped oscillations in I c as a function of <1f and ex- 
change energy in the clean 28 and dirty limits. 29,30 These 
changes in the critical current have been experimentally 
confirmed 10-18 . Of particular interest is Rcf. 16, which 
demonstrates the robustness of <-> n transitions by pro- 
viding evidence of switching in samples with interfaces 
that were not atomically smooth. Indeed, despite de- 
viations as large as 0.6nm over 10% of a sample, clear 
evidence of switching was found. Near T c , and in the 
diffusive limit, the theory was later extended to include 
arbitrary interface transparency. 30 Measurements of the 
superconducting phase 13 have corroborated the n state 
in SFS junctions, and the predicted oscillations in several 
thermodynamic quantities have in many cases been found 
experimentally. Direct evidence of DOS oscillations was 
reported in a tunneling spectroscopy experiment, 31 but 
not observed 32,33 in other cases. Such studies give us the 
valuable insight that the oscillations are correlated with 
7r <-» transitions. In this work, we show that there is 
indeed an intimate relation between the oscillations in T c 
as a function of relevant parameters and the transitions 
from the n to the state and we find good quantitative 
agreement with experimental data. 

Since the possibility of having a particular junction 
configuration depends fundamentally on the intricate 
properties of the pair amplitude, the complicated and 
demanding task of calculating the pair potential, A(r), 
rigorously and self-consistently becomes absolutely nec- 
essary, particularly as the inhomogencities occur on a 
microscopic scale. The first step in the self-consistency 
process often involves an assumed simple piecewise con- 
stant form for A(r), which is then iterated through the 
relevant equations until convergence is achieved. It is 
not justified to bypass the technical difficulties associ- 
ated with self-consistency and to use only an assumed 
form for the pair potential. The final calculated A(r) of- 
ten deviates significantly, even in overall symmetry, from 
the assumed form: the self-consistent A(r) has a compli- 
cated spatial behavior that can lead to stable states mix- 
ing and 7r junction configurations. 34 A self-consistently 
calculated pair potential minimizes, at least locally, the 
free energy of the system. To determine whether the 



calculated state is merely a local minimum of the free 
energy or the global one, the free energies from all possi- 
ble self-consistent 0- and 7r-junction configurations must 
be compared with high precision. Recently developed 
numerical algorithms 34,35 overcome the difficulties that 
arise in computing the small difference between much 
larger quantities and enable accurate computation of the 
differences in the values of the condensation free energy 
of different minima. 

For clean SFS junctions, a relevant set of basic pa- 
rameters to consider includes dp, the exchange energy 
ho and T. As these parameters vary, the or 7r-state 
free energies may cross at certain points in parameter 
space, yielding phase transitions. It has been shown 34 
that at T = 0, transitions occur when varying ho, dp 
and also the mismatch parameter A, defined as the ratio 
of Fermi energies in the F and S regions. This mismatch 
can induce a transition because at A w 1, when the Fermi 
wavevectors match, the layers couple more strongly, while 
at small A the coupling is effectively weaker. If the tem- 
perature varies it is also possible to have a first order 
transition between and tt junction states, as recently 
shown in both the clean, 35 and dirty 36 limits, and also 
predicted for short-period F/S superlattices. 26 The tem- 
perature has been shown to have a pronounced effect on 
the pair amplitude in the F region of F/S structures, 37 
strongly diminishing its magnitude while maintaining its 
characteristic period of oscillation as T increases. This 
translates into weaker coupling between adjacent S lay- 
ers. If the magnet width is such that the junction is 
near a <-» tt transition point at T = 0, increasing the 
temperature can result in the critical current of the junc- 
tion having a non-monotonic temperature dependence. 25 
It has been argued 38 that the transition is discontinuous 
in uniform samples but rounded off in samples of vari- 
able thickness. 39 However, a transition can be observed 16 
in just a portion of samples with nonuniform thickness. 
These results indicate that the temperature can be used 
to switch between a and a tt state configuration. It is 
possible to locate regions of parameter space that give the 
desired transitions using the T = results as guides, how- 
ever the task is still significantly demanding. Such tem- 
perature transitions were found to occur in one-junction 
and 3-junction systems for moderate values of A. 35 Thus, 
a 1-junction system was found to have a — > tt first order 
transition as T was lowered, and a tttttt — * ttOtt transi- 
tion was found for a 3-junction system. 35 In each case, the 
free energies of a stable and a metastablc state crossed at 
the transition temperature with differing derivatives, and 
therefore entropies. The existence of metastable states 
and an entropy discontinuity arc hallmarks of first-order 
phase transitions. Moreover, the reported latent heats 
were reported to be within available experimental reso- 
lution. It is therefore desirable to systematically study 
the coexistence of metastable states and the nature of the 
transition in SFS and higher-order multilayer structures. 

The main objective of this paper therefore, is to map 
out the regions of parameter space in which the different 
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junction states are stable, and to trace the locations of 
the phase transitions in systems with SFS junctions. An 
extensive sweep of the geometric and materials parame- 
ters including d F , A, as well as T, is performed. To start 
with, it is important to know which d F and A ranges al- 
low more than one self-consistent state at T = 0. One 
can then check if a metastable state at low temperature 
becomes the equilibrium state at higher T. By using 
this procedure we obtain a complete phase diagram of 
an SFS junction within the relevant region of (T, A, dp) 
space. To accomplish this, we use a method that can 
accommodate arbitrary values of the above parameters, 
without recourse to approximations. As discussed above, 
all calculations involving the pair potential must be per- 
formed using fully self-consistent algorithms, starting 
from the microscopic equations (Bogoliubov-deGcnncs 
(BdG)). The need for a fully microscopic theory arises be- 
cause the characteristic period of the pair potential oscil- 
lations approaches the atomic scale. For the nanoscale in- 
tcrlaycr widths considered here, geometrical oscillations 
decisively influence the final results. 

We present in Sec. II the microscopic equations and 
the associated notation relevant for systems containing 
SFS junctions. We review the numerical procedures in- 
volved in calculating the self-consistent pair potential 
and quasiparticle spectra, and the method used to cal- 
culate the primary thermodynamic quantity, the conden- 
sation free energy, AT(T), from the self-consistent spec- 
trum and pair potential. We also outline a semi-analytic 
method to calculate T c through the linearized BdG equa- 
tions. In Sec. Ill, we show that first order transitions 
with measurable latent heat can occur between states 
containing different numbers of and tt junctions as the 
temperature changes. For SFS junctions the transitions 
we find are from the tt to the state as T increases, as 
found in experiment, 16 and occur predominantly in re- 
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Here, H = p 2 z /2m — Ep {z) + e± is a single-particle Hamil- 
tonian where p z /2m + e±_ is the kinetic energy term. The 
continuous variable ej_ is decoupled from the z direction 
but of course it affects the eigenvalues e n . We describe 
the magnetism by an exchange field h(z) which takes the 
value ho in the F material and vanishes in S. Within the 
superconducting layers, E F (z) is equal to Ep$, the Fermi 
energy of the S layers measured from the bottom of the 
band, while in the ferromagnet we have E F (z) = Epm 
so that in the F regions the up and down band widths 
are Ep-\ = Epm + h and Epi = E F m — ho respectively. 
In the seven layer case we assume parallel orientation of 
the magnetization in all F layers. One should not as- 



gions where T c is low. Using the T c calculated from the 
linearized theory and the <-» tt phase transitions, we 
obtain the full phase diagram in an extended region of 
parameter space spanned by T, dp, and A. We com- 
pare our calculated oscillations in T c as a function of dp 
with reported Nb/Co experimental data 20 and find good 
agreement. 



II. METHODS 

The systems that we study consist of slabs of clean 
superconductor (S) material separated by ferromagnetic 
(F) layers. We will emphasize trilayers consisting of one 
SFS junction and, as a sample of what can generally oc- 
cur in multilayers, present also results for seven layer 
systems consisting of three SFS junctions. The thickness 
of the S layers in an SFS junction is denoted by ds, and 
that of the F layers by dp. The seven layer system con- 
sists of three SFS junctions stacked together, so that the 
thickness of the two inner S layers is 2ds- We assume that 
the layers are semi-infinite in the directions perpendicu- 
lar to the interfaces (the x — y directions) and that the 
interfaces are sharp. The spatial inhomogeneity is con- 
fined to the z direction, allowing us to model the system 
as quasi one dimensional. We assume parabolic bands, 
thus in the transverse direction e± = k\ /2m, where k± 
is the wavevector in the transverse direction and e± is 
the energy corresponding to the x — y variables. 

We use the microscopic Bogoliubov-deGennes 40 equa- 
tions to study this inhomogencous system. Given a pair 
potential (order parameter) A(z) that is to be deter- 
mined self consistently, the spin- up quasi-particle (u^z)) 
and spin-down quasi-hole {v+(z)) amplitudes obey the 
BdG equations in the following form: 



(2.1) 

I 

sume that E FM = E FS and we therefore introduce the 
dimensionless Fermi wavevector mismatch parameter A 
by Epm = AEps- Usually, one has A < 1. It is also con- 
venient to introduce the dimensionless magnetic strength 
variable / by ho = EfmI- The 1 = 1 limit corresponds 
to the "half-metallic" case. We neglect interfacial scat- 
tering. The amplitudes u^z) and v\{z) can be written 
down from symmetry relations. 40 

The required self-consistency condition for the pair po- 
tential A(z) is: 

n 

(2.2) 



H - h(z) 
A(z) 



A(z) 
<H + h(z)) 



u\{z) 
vj,{z) 
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where here and below the prime indicates a summation 
over states for which |e„| < cup, where urj is the usual 
cutoff "Debye" energy and it is understood that the index 
n includes k± as well as the longitudinal variables. The 
BCS coupling g(z) is taken to be a constant g in the 
superconductor and zero in the fcrromagnct. 

A. Self-Consistent Solutions 

Equations (2.1) and (2.2) comprise a non-linear set of 
equations. Exact solutions to this set must be computed 
in a self-consistent manner. We follow the procedure 6 ' 7 ' 35 
used in previous work; we omit the repetition of the tech- 
nical details. We begin with an assumed form for A(z), 
either from a prior calculation with similar parameters or 
an a priori guess (usually a stepwise function) , and then 
numerically solve Eq. (2.1) for every value of in the 
appropriate range to compute u^z), v^(z), and e„. An 
expansion of all quantities in terms of sine waves is used 7 
to carry out the solution. The required matrix elements 
are given, for our geometries, in Ref. 6. This resulting en- 
ergy spectrum and quasiparticle amplitudes are used in 
Eq. (2.2) to compute a new A(z). We then feed this new 
A(z) back into Eqs. (2.1) and repeat this process until 



T(T) = -2T^'ln [2 

n 



where d is the total thickness of the system in the z- 
direction. In this expression, only the energy eigenvalues 
appear explicitly, the cigcnfunctions appearing only in- 
directly through the self-consistent A(z). It is equivalent 
to several other expressions found in the literature which 
contain the quasi-particle amplitudes explicitly, but it is 
computationally much more convenient. 

The condensation free energy is defined as AT(T) = 
Ts—3~n, where Ts is the free energy of the superconduct- 
ing state and that of the non-superconducting sys- 
tem. We compute by setting A = in equations (2.1) 
and (2.3). Calculating AT(T) is a significant numer- 
ical challenge: recall that in a bulk superconductor 42 
AF(0) = -(1/2)7V(0)A 2 , where N(0) is the usual den- 
sity of states and A is the order parameter for the bulk 
superconductor at T = 0, which is several orders of mag- 
nitude smaller than Tn oc N^ujq. Hence, to obtain AT 
we must subtract two numerically obtained large quanti- 
ties in order to extract a difference several orders of mag- 
nitude smaller than the terms subtracted. Furthermore, 
as we shall see, the difference in condensation free ener- 
gies of competing self-consistent states (when they occur) 
is only a small fraction of the condensation free energy 
of each of them. To obtain sufficiently accurate values 



the fractional difference between the average of succes- 
sive solutions for A(z) is less than a threshold value that 
we take to be 10~ 5 . Solutions obtained in this way are 
exact up to the chosen numerical precision. 

The self-consistent solution for a trilayer SFS junction 
can be of the it or the type, with the pair potential 
either changing or not changing sign across the F layer, 
respectively. More complicated situations can occur in 
multilayers: for a three junction system one can en- 
counter four symmetric states (000, 07r0, irOir, nmr, with 
each symbol corresponding to the state of each junction). 
When, for a given temperature and set of geometrical and 
material parameters such as /, (If, and A, several 6 ' 7 ' 35 
different self-consistent solutions, that is, local minima 
in the free energy, exist, the stable state must be deter- 
mined by comparing the condensation free energies of the 
competing self-consistent states. As discussed in Sec. I, 
when the equilibrium state changes as a function of tem- 
perature 35 a first order phase transition can occur, with 
a corresponding latent heat. One of the chief goals of this 
paper is to study an extended region of parameter space, 
locating where such transitions exist and then mapping 
out the corresponding phase diagram. 

To evaluate the free energy, T, of the self-consistent 
states we use the formula from Ref. 41: 
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of AT requires therefore a very high degree of precision 
in calculating Ts and Tn so that we can distinguish the 
relatively small differences between competing states to 
locate phase transitions. This situation is made more 
challenging by the need to calculate derivatives of AT to 
obtain thermodynamic functions and latent heats. 

B. Calculation of T c : Linearized Solution 

While the transition temperature T c from the non- 
superconducting to the superconducting state can be nu- 
merically calculated as the temperature at which AT 
vanishes, it is much easier to evaluate T c by treating 
A(z) as a small parameter and linearizing the equations. 
In this way the calculation is nearly entirely analytic. 
The amplitudes are written as u^z) — u 6 (z) + u' n (z) 
and v^(z) — v®(z) + v' n (z) (we have dropped the spin 
indices for simplicity). The u^(z) and v®(z) terms arc 
computed from the zeroth order equation, which is ob- 
tained by setting A(z) = in Eq. (2.1). The form of the 
zeroth order equation implies that u°(z) and v®(z) are 
completely decoupled and have distinct energy spectra, 
denoted by and ejj respectively. Proceeding to calcu- 
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late the lowest order corrections, we incorporate quasi- 
particle coupling through the pair potential matrix. One 
can then obtain u' n (z) and v' n (z) from textbook pertur- 
bation formulas. The intermediate sums are in principle 
over the entire zeroth order spectrum, but as a practical 
matter it is enough to include in these sums energies 
and within a few ud of the Fermi level. 

We then expand the quasiparticle amplitudes and their 
first order corrections in a sine wave basis 4> q (z), e.g. 

U n{ Z ) = J2q Uqn<j) q (z) and v' n (z) = J2q V qn ( l ) q( Z ), WnCrC 

4> q (z) = \j2jd sin(fc 9 z), with k q = qn/d. The range 
of the sums over k q is formally infinite, but again it is 
only necessary to sum up to a wavenumber fcjy with an 



associated energy a few Lop, from Ep. Inserting these ex- 
pansions into Eq. (2.2) gives the lowest order correction 
to A(z), which we then expand in the 4> q (z) basis. Upon 
taking into account the orthogonality of the basis func- 
tions, the expanded Eq. (2.2) is then transformed into 
the matrix equation 

A ; = J ik^k, (2.4) 
k 

where the are the expansion coefficients of A(z) in 
terms of 4>k{z). One finds after straightforward algebra: 



Jik 




pql' 



T N v" U° K h 



U° K 

pn qm '■ 



pql' 



m pq 



tanh ( ^ ] + 



tanh 




(2.5) 



Here we have used gKijk = J g{z)4>i{z)cj)j{z)(j)k(z)dz. 
The integral over e± reflects the dependence of the zeroth 
order quasiparticle amplitudes and energies on e±, and 
the sum over n is here only over longitudinal quantum 
numbers with the prime denoting the limitation indicated 
below Eq. (2.2) on the energies e£ and 

The transition temperature can then be found, in anal- 
ogy with standard procedures, 43,44 by treating Eq. (2.4) 
as an eigenvalue equation for the matrix J^. At the 
transition temperature T c the largest eigenvalue is unity, 
while if T > T c all eigenvalues are less than unity. Un- 
like the free energy method described in subsection II A, 
this procedure does not require an iterative process and 
only the last step (finding the eigenvalue) must in prac- 
tice be performed numerically. Therefore this method is 
much more efficient, and it also provides a check on the 
numerics of our free energy. 



Similarly, we have the average number density for each 
spin subband, 

<M*)> - E {KM! 2 + K«] 2 1 1 - • 

n 

(2.7) 

0. 

(2.8) 



This leads to the dimcnsionlcss magnetization, M(z), 
(n T (z)) - (»»;(*)) 



M(z) 



(ni( Z )) + ( ni (z))' 



which reduces to M(z) = 

[(1 + If' 2 - (1 - If/ 2 ] I [(1 + If I 2 + (1 - If/ 2 ] 
for bulk F material, within our assumptions. This 
expression is in numerical agreement with our results 
from Eq. (2.8) in sufficiently thick F layers. 



III. RESULTS 



C. Other quantities 

From the self-consistent amplitudes and an energy 
spectrum we can also calculate other quantities of in- 
terest such as the density of states (DOS) and the mag- 
netization. The local density of states is 

N(z, e) = -EE [K W] 2 - e «) + K W] 2 /'( £ + e ") 

a n 

(2.6) 

where a denotes spin and /'(e) is the first derivative of 
the Fermi function. One can also omit the sum over a 
and obtain the spin dependent DOS. 



In this section we present and discuss our results. As 
explained above, trilayers consisting of one SFS junction 
will be emphasized but results for seven layer systems 
comprised of three stacked SFS junctions will also be 
given in order to show that the single junction results 
can indeed be generalized to multilayer samples. We will 
first discuss the results for the thermodynamics and the 
phase transitions that ensue. This will include a detailed 
discussion of the phase diagram for the SFS trilayer in 
the most interesting region of the three dimensional space 
spanned by T, A and the F layer thickness. A discussion 
of the properties of the transition temperature T c as a 
function of dp and a comparison with experiment follow. 
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We will also discuss other quantities of interest, such as 
the pair amplitude, the density of states, and the mag- 
netization. 



A. Parameters and Units 

The results presented below will be given in terms of 
convenient dimensionless quantities. We measure all the 
lengths in units of kpg, the Fermi wavevector in S. We fix 
Ds = kpsds = 100. We have taken the BCS coherence 
length £ equal to ds- For ds of order of or larger than £ , 
results are only weakly dependent on ds, hence our re- 
sults are applicable to a very wide range of values of this 
variable, provided ds is not too small. The dimensionless 
thickness Dp = kpsdp of the ferromagnetic layers will 
be varied over the range of interest, which corresponds to 
relatively small values, since at large ones the F/S prox- 
imity effects are negligible. Similarly, we introduce the 
notation Z = kpsz. The magnetic strength parameter is 
taken to be / = 0.2 unless otherwise noted. The effects 
of varying / are physically similar to those of varying Dp 
since the pair amplitude oscillations in F are governed 4 
by the difference — ki)dp between Fermi wavevectors 
in the spin bands in F. The Fermi wavevector mismatch 
parameter, A, to which results are quite sensitive, is var- 
ied over the experimentally relevant range 0.1 < A < 1. 
The temperature is measured in units of T°, the criti- 
cal temperature of a bulk sample of the material S. We 
choose ujd/Eps = 0.02; the dimensionless quantities cal- 
culated are not sensitive to this choice. Condensation 
free energies will be given in units of -/V(0)Aq, twice the 
absolute value of the condensation free energy of a bulk 
superconducting sample of the same total thickness at 
T = 0. A dimensionless measure of the latent heats will 
be given by dividing the corresponding entropy discon- 
tinuities by C n (T®), the specific heat of a sample of the 
same overall thickness but consisting exclusively of the S 
material in its normal metal state at T°. 

We performed several checks of our numerical meth- 
ods. We verified that the temperature at which the 
self-consistent condensation free energy goes to zero is 
in each case the same as the transition temperature ob- 
tained from the linearized solution to the BdG equations. 
For a sample with ds 3> £o and dp — 0, we quantitatively 
recovered the well established results 42 for the thermody- 
namics, including the second order phase transition at T° 
and the associated specific heat discontinuity. Further- 
more, the spatially averaged DOS computed numerically 
for this system shows a well defined gap at energies within 
A of Ep and the characteristic divergence at Ep ± A . 
This test is very severe since our numerical method must 
necessarily be more accurate for smaller systems, where 
fewer variables are required. Thus the ability of our nu- 
merical procedures to handle the relatively large systems 
(over six superconducting correlation lengths thick) con- 
sidered here is verified. The low temperature limit was 
extensively checked in Ref. 6, and it was also previously 



verified 7 that our methods give the correct thickness de- 
pendence of A(z) for a superconducting slab as found in 
the literature 45 . 



B. Free energy 

The basic quantity that determines the phase transi- 
tions and the entire thermodynamics is the condensation 
free energy, AJ-(T). We will therefore begin our expo- 
sition by describing some of our results for AjF(T) at a 
few parameter values representative of the regions where 
the phase transition behavior is richest. 

We begin with Fig. 1 which shows, for an SFS trilaycr, 
the self-consistent condensation free energy AjF(T) plot- 
ted versus reduced temperature T/T®. Data points were 
obtained at T/T c ° = 0.01 intervals. The values of Dp 
and A for which results are shown were chosen to be 
such that, at T = 0, self-consistent solutions of both the 
and the it states are found to exist. In each panel, the 
free energies of the two competing states are shown. The 
thermodynamically stable state is of course the one with 
the lower free energy. The slope of AT(T) approaches 
zero as T — > 0, which indicates that the calculated en- 
tropy vanishes at T — as required by the third law 
of thermodynamics. The slope also approaches zero as 
Af{T) — > 0, indicating that the transition to the nor- 
mal state is of second order, without latent heat. The 
temperature at which this second order phase transition 
occurs, T c , is the temperature at which the lower AT(T) 
vanishes. The T c found this way agrees with the indepen- 
dently calculated T c using the linearized BdG equations. 
The inherent finite-size and proximity effects cause T c to 
be considerably smaller than T° in all cases. 

Two outstanding features of the results shown in this 
figure are the existence of a metastable state at all tem- 
peratures from zero up to T c and a first order phase tran- 
sition at an intermediate temperature: the free energies 
of the two competing states cross at T — T x . In all cases 
shown the tt state is stable below T x and the state is 
stable above. The position of T x is marked by the vertical 
arrow in each panel. The value of T x changes smoothly 
as Dp or A are changed. One can see by comparing top 
and bottom panels how T x changes with Dp at constant 
A. 

In Fig. 1 one can see the difference in slope between 
the stable state and the metastable state at T x , partic- 
ularly apparent in the right panels. The existence of a 
metastable state and the discontinuity in the slope of the 
free energy of the stable state (i.e., the entropy) at T x 
indicate the existence of a first order phase transition 
with an associated latent heat. That any such transi- 
tion should be first order can also be expected from the 
change of symmetry of the pair amplitude. For a range 
of parameter values including those shown in this figure, 
the phase transition behavior is exceptionally rich. In 
many other regions of A and Dp parameter space the 
behavior is simpler: in some there is only one self con- 
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FIG. 1: (Color online) Results for the normalized (see text) condensation free energy AT(T) vs. temperature for a 3 layer 
SFS junction. The different curves are labeled in the legends. In all cases shown, upon increasing T a ir to transition 
occurs at temperature T x , indicated by the arrows. The top left panel shows results for A = 0.550 and D F = 7.0, resulting 
in T x /T° = 0.07. Bottom left: A = 0.550 and D F = 7.1, with T x /T° = 0.16. Top right: A = 0.650 and D F = 5.8, with 
T x /T° = 0.13. Bottom right: A = 0.650 and D F = 5.9, resulting in T x /T° = 0.23. 



sistent solution to the BdG equations at T = 0, while 
for other ranges of A and Dp a metastable state is found 
at low temperature but it never becomes the stable state 
as T increases. It is only in some regions of parameter 
space that «-» 7r transitions occur as a function of T. 
This question will be discussed in more detail below. 

Examples of similar results for the 7 layer case are 
shown in Fig. 2. These are all at A = 0.55, for several 
values of Dp. In all cases shown at least two of the four 
possible metastable states mentioned above exist over the 
entire temperature range. The states shown in each panel 
are the two lowest in free energy. In some cases additional 
states exist but with higher free energy throughout: any 
such states are omitted from the plots. The three panels 
illustrate three different types of phase transitions as T 
increases: nOn — > 7r7T7r (one junction flipping — > n), 
OitO — > 7r07T (three flips), and 07r0 — > 000 (one flip, ir — » 
0). Each of these persists over a range of Dp . The results 
show all of the same qualitative features as the three layer 
case: the slope of AJ-(T) approaches zero as T — ► and 
as AT(T) — > 0. There is again a change in the slope of 
the stable state at T x which shows that the transition is 
also of first order in multi-junction cases. We will show 
that the latent heats are of the same magnitude as or 



larger than in the 3 layer case. There is an important 
quantitative difference: T x varies more slowly with Dp or 
A and therefore the range of parameter values for which 
such transitions are found is wider. One can expect then, 
that in higher order multilayers these phenomena will be 
even more general. 



C. Thermodynamic functions 

From the free energy one can obtain the entire thermo- 
dynamics. Figure 3 shows some of the thermodynamic 
functions that can be obtained from the results shown in 
Fig. 1. Results are shown for two quantities: the dimen- 
sionless condensation entropy S(T), defined as the neg- 
ative derivative of AT(T) with respect to the reduced 
temperature T/T®, and the dimensionless condensation 
energy U(T) defined as U(T) = AT{T) + (T/T C °)S(T). 
Results are shown for both the stable and the metastable 
states as a function of reduced temperature. One sees 
that S — > smoothly as T — > and T — * T c in each case, 
which is an important check on the computation. The 
condensation entropy, energy, and free energy all vanish 
at T c . In each panel a bold vertical line indicates T x . 
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FIG. 2: (Color online) Results for A.F(T) vs. reduced tem- 
perature (as in Fig. 1) for a 7 layer system. T x is indicated by 
the arrows. The different curves are labeled in the legends. 
The top panel shows a 7r07r — > nnn transition for A = 0.55 
and D F = 9.1. The transition occurs at T x /T° = 0.37. The 
middle panel shows a 07r0 — » 7r07r transition for A = 0.55 
and D F = 7.9, at T x /T° = 0.33. The bottom panel shows 
a more pronounced 07r0 — > 000 transition for A = 0.55 and 
D F = 4.75, with T x /T° = 0.27. 



The free energy crossings correspond neither to crossings 
in S(T) nor to crossings in U(T). The former follows 
from the phase transitions being of first order with an 
associated discontinuity in the entropy. Both the energy 
and the entropy, therefore, play important roles in the 
phase transition. The specific heat is not shown but can 
be calculated by taking a further derivative. 
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FIG. 3: (Color online) Thermodynamic functions of an SFS 
trilayer. The top and bottom panels show the condensation 
energy and entropy (in dimensionless form, see text) for the 
two sets of parameter values used, respectively, in the left 
and right bottom panels of Fig. 1. The meaning of different 
curves is indicated in the legend. The location of the first 
order transition is marked by the bold vertical line. 



Examples of the thermodynamic functions for the 7 
layer system are shown in Fig. 4. The 3 junction case 
is again qualitatively much like the one junction case. 
The entropy, energy, and free energy all go to zero in the 
appropriate limits. There is no crossing in energy or en- 
tropy at T x , indicating the important interplay between 
energy and entropy. The discontinuity in the entropy at 
T x also reflects a latent heat, comparable to or larger 
than that in the one junction case. 

The behavior of the Cooper pair amplitude at the first 
order transition is illuminating. Figure 5 shows, for the 
SFS trilayer at T = T x , the pair amplitude (defined in 
the usual way as the average of spin up and down cre- 
ation operators) normalized to Ao/g, its value in bulk 
S material at T = 0. Results are given versus dimen- 
sionless position Z. The F region is in the middle, set 
off by vertical dotted lines, and only small portions of 
the S regions are shown. The two cases shown corre- 
spond to the two bottom panels in Fig. 1. The absolute 
value of the pair amplitude is discontinuous at T x : in 
both plots it is slightly larger for the state. We recall 
that for a bulk superconductor at T — 0, the free energy 
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FIG. 4: (Color online) Thermodynamic functions of a 7 layer 
system. The panels show the condensation energy and en- 
tropy (as in Fig. 3) for the parameter sets used in Fig. 2. The 
different curves are labeled in the legend, extending the no- 
tation introduced in Fig. 3. The location of T x is marked by 
the bold vertical line. 



is proportional to the average value of the squared pair 
potential, 42 and this is also 6 approximately the case at 
T = for SFS layered systems when d,F <C ds <C Co- 
Even in the bulk system, however, such a relationship is 
not valid 46 at finite temperature. It is therefore unrea- 
sonable to expect this proportionality in the layered case 
and indeed it does not occur: the pair amplitudes do not 
change continuously at T x . We conclude that the phase 
transitions are not driven by the pair amplitude. 

The pair amplitude for the three junction system dis- 



FIG. 5: (Color online) The normalized (see text) pair ampli- 
tude for the and n states of the SFS trilayer, at the crossing 
point T x , as a function of position Z = kFSZ. Only the mid- 
dle portion of the sample is shown. The F layer is delimited 
by the vertical dotted lines. Results are presented (top and 
bottom panels) for the two sets of parameter values used, re- 
spectively, in the left and right bottom panels of Fig. 1. 



plays properties that are very similar to those of a sin- 
gle junction. A representative example, corresponding to 
A = 0.55 and Dp = 4.75, is shown in Fig. 6. The abso- 
lute value of the amplitude is again discontinuous at T x . 
It is very important that the 7 layer and the 3 layer sys- 
tems have qualitatively similar properties, as this shows 
that the phenomena we discuss are very general. At the 
same time, in the 7 layer case there is a greater number 
of possible transitions and the regions of parameter space 
in which they occur as a function of T are wider, indi- 
cating that such phenomena can be more readily found 
in more complicated systems. We can make qualitative 
predictions for the 7 layer system based on our quanti- 
tative (but computationally less demanding) calculations 
for the 3 layer system. Thus, the properties of a single 
junction system can be generalized to systems with many 
junctions. 



D. Latent heats 

The signature of a first order phase transition is its 
latent heat. In Fig. 7 we show results for the dimen- 
sionless latent heat i, defined as the difference between 
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FIG. 6: (Color online) The dimensionless (see text) pair am- 
plitude for the 07r0 and 000 states of the 7 layer system for 
A = 0.55 and Df = 4.75 at the crossing point T x as a function 
of position Z. This is the parameter set used in the bottom 
panel of Fig. 2. 



a factor of two, than the upper end of the scale. 

The latent heats vanish as T x approaches or T c , con- 
sistent with the smaller condensation entropy of each 
state in those limits. However, whenever T x does not 
approach these limits the latent heat can exceed 1% of 
C„(T C °) for one junction, and even more for the three 
junction system. Since we give L in units of C n , which is 
an extensive property, it should be easier to observe these 
latent heats in larger systems. A value of L w 0.01 would 
correspond to picojoules in actual samples of relatively 
small size. 10 Such latent heats can be readily observed via 
standard techniques used to measure specific and latent 
heats in films. 47 Even smaller specific heats can be mea- 
sured using multiple samples: attojoule level results have 
been reported 48 in electronic systems. We see therefore 
that whenever a first order transition occurs, the associ- 
ated latent heat is observable. 
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FIG. 7: Latent heats. L is the entropy discontinuity in units 
of C n (Tc) (see text). It is plotted against the reduced temper- 
ature of the first order phase transition. The symbols joined 
by lines are for an SFS trilayer: the value of T x is changed 
along the horizontal axis by varying Df, and from curve to 
curve by varying A (see legend). The triangles are for the 7 
layer system cases shown in Fig. 2. The vertical arrow at- 
tached to the topmost triangle indicates that it corresponds 
to a value L = 0.029 (off the scale). 



the entropy of the stable states just above and just be- 
low T x divided by C n (T°), the specific heat of a normal 
bulk sample of S material at T — T®. This is appro- 
priate because C n (T) is equal to the entropy in the free 
electron model. Results are plotted as a function of T x . 
Most of the results shown are for a single junction: in 
that case the crossing temperature is varied by changing 
Dp for several different values of A, as indicated by the 
symbols connected by straight segments. The three data 
points indicated by the isolated triangles correspond to 
the three transitions shown in Fig. 2 for the 7 layer sys- 
tem. One of them corresponds to a value larger, by over 



E. Phase diagram 

We have seen that in an SFS trilayer there are two 
kinds of phase transitions. First, there are second order 
phase transitions from the normal state to a supercon- 
ducting state of either the or the 7r kind. There are also, 
at certain ranges of the relevant parameters, first order 
transitions between the and 7r superconducting states. 
As a practical matter, observability of the latter transi- 
tions through thermodynamic measurements requires an 
appreciable difference in condensation energies between 
the two states. This difference is an oscillatory function 
of Dp at constant A and / (see e.g. figure 3 of Ref. 6) 
with the oscillations becoming damped at large Dp, since 
then, at any / > the proximity effects are reduced and 
the and ir states are degenerate. Hence the most im- 
portant regions theoretically and experimentally are at 
relatively small values oi Dp. As to A, the entire region 
A < 1 is relevant. 

Therefore, we have mapped out the entire phase dia- 
gram of an SFS trilayer in this most relevant region of 
(T, A, Dp) space in Fig. 8. As explained above, vary- 
ing / is equivalent to varying Dp, so we use Dp as the 
more experimentally relevant parameter. We show two 
views of the phase diagram to aid in the visualization of 
this three dimensional figure. There are three regions in 
this diagram, each representing one of the three possible 
states: state, n state, and normal (not superconduct- 
ing) state. The crossings T x are calculated from the free 
energies, and T c through the linearization method. 

The top sheet shows the superconductor/normal metal 
transition. As Dp — > 0, T c /T® approaches unity for all A. 
At small A the sheet also flattens, since then the Fermi 
level of the ferromagnet is small compared to Eps so that 
there is little interaction between the Cooper pairs and 
the ferromagnet. The finite temperature T x transitions 
between and it regions are located at the sheet or "wall" 
that goes from the T = plane to the T c sheet, separating 
the from the w state regions. This wall is of course not 
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FIG. 8: (Color online) The (A, Df,T) phase diagram for the 3 layer system. The two panels show different views of the same 
plot. There are three regions: in those labeled and n, the and n states are, respectively, the equilibrium state, while normal 
metal indicates where the sample is nonsuperconducting. The top surface separates non-superconducting and superconducting 
regions. The fairly vertical sheet marks the temperature transitions between and ir states. The intersection of the — w and 
the T c boundaries is marked by a dotted line. The portion of the T = plane marked by x symbols is the range of (A, Df) for 
which only the n state exists for all T: there is no metastable state of the type. Likewise, in the region marked by + symbols 
only the state exists. In the portion left blank, solutions of both kinds are possible. 



completely vertical: its deviation from verticality is what 
causes first order phase transitions as a function of T. On 
the smaller Dp end, this wall ends because one of the two 
states becomes unstable: a region of parameter space is 
entered where only one self consistent solution exists at 
any temperature. Coincident with this, as one can see 
more clearly in the right panel, T c is sharply reduced: 
in other words, the condensation energies of both states 
rise towards zero, with one actually vanishing. Near this 
region T c has always a sharp dip. As one proceeds to- 
wards the opposite end of the wall, at larger values of 
Dp, T c increases and the wall becomes steeper, until it 
eventually becomes vertical. Beyond that, no transition 
occurs as a function of T: the stable state is the same 
at all temperatures. Beyond the portion shown, there- 
fore, the wall would become completely vertical and it is 
not depicted because it would obscure the diagram. It is 
sufficient to show its behavior in the T = plane. The 
crossings at T = are not thermodynamic phase transi- 
tions, they merely indicate a change in the stable state 
as various sample parameters are changed. 

Exploration of T c and AJF(O) for larger values of Dp 
at several values of A indicates the existence of other 0- 
7T boundaries at larger values of Dp. Thus, one could 
extend the phase diagram in that direction, but as previ- 
ously seen 6 and discussed above, these additional regions 
are qualitatively similar to the one shown here in detail, 
and quantitatively less interesting. 

Computing a complete three dimensional phase dia- 
gram such as the one in Fig. 8 for a 7 layer system 
would be very expensive in computational resources and 
is not necessary. We have already seen in connection 
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FIG. 9: The T = plane of the phase diagram for the 7 layer 
system. The regions in which each of the four possible sym- 
metric states is the stable one are indicated by the symbols in 
the legend. There are also metastable states at most values 
of A and Df- 



with Fig. 2 that first order transitions not only occur but 
are more abundant in such systems. Further evidence is 
in Fig. 9 where we show the zero temperature plane of 
the 7 layer phase diagram. A different symbol marks re- 
gions where each of the four possible symmetric states is 
the stable one at very low T. The many boundaries be- 
tween the various states and the presence in many areas 
of metastable states (not marked) reflect that there may 
be many first order phase transitions in the 7 layer sys- 
tem. We can thus infer that even larger structures will 
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FIG. 10: T c vs. D F (left scale, * symbols) for A = 0.70. Also 
shown (right scale) are the energies of the and tt states at 
T = (+ and x signs respectively). Note the correlation 
between the changes in the stable state at T = and the dips 
in T c . 
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FIG. 11: Calculated values of T c , in good agreement with 
experimental 20 data for a Nb/Co system. The experimental 
points shown are the average of the two series reported in 
Ref. 20. 



have rather intricate and rich phase diagram we found 
that the range of parameter space over which each type 
of transition persists is much broader than in the 3 layer 
case. For example, for a fixed A = 0.55, we observed 
phase transitions (see Fig. 2) between 000 and 07r0 states 
for 4.4 < Dp < 5.0, transitions between 07r0 and ttOtt 
for 7.75 < Dp < 8.0 and transitions between irOir and 
tttttt for 9.0 < D F < 9.5. For a fixed Dp = 10, the 
ttOtt <-> tttttt transition exists 35 for 0.35 < A < 0.50. We 
did not search for other transitions at Dp = 10. 

By taking a slice of the phase diagram in Fig. 8 at 
fixed A, one can discern regular, damped oscillations of 
T c with Dp. In Fig. 10 we show T c for A = 0.70 over 
an extended range of Dp. It is clear that as Dp is in- 
creased the amplitude of the oscillations decreases. This 
is in good agreement with experiment, as we shall see in 
detail below. In addition to T c , this figure shows AJ r (0) 
for the and tt states. In a bulk superconductor, the 
ratio of this dimensionless quantity to the reduced tran- 
sition temperature is —0.5, which is confirmed here by 
our result for the state at Dp = 0. The tt state is 
unstable, for obvious reasons, in the Dp — ► limit. At 
finite values of Dp this relationship between normalized 
condensation energy and reduced transition temperature 
is not strictly obeyed, but there is a qualitative correla- 
tion: increases in the absolute value of A^ r (0) correspond 
to increases in T c . The values of Dp at which the stable 
state switches between and tt correspond to the sharp 
dips in T c in all cases. This has also been seen in connec- 
tion with Fig. 8 and it indicates that the structure and 
shape of the oscillations in T c are strongly correlated with 
the low temperature state. The free energy data plotted 
have gaps, notably for the tt state near Dp < 4 and for 
the state at 8 < Dp < 17. These values of Dp delimit 
regions in which the self-consistent calculation resulted 
in only one state. The free energy of the vanishing state 
goes continuously to zero at those boundaries. The pair 
amplitude is found to go also smoothly to zero. 



The low temperature crossings at the many different 
values of Dp suggest the location of more <-> tt phase 
transitions. This is in agreement with the direct observa- 
tions reported in Ref. 24. Another corroboration of this 
claim comes from Ref. 16, in which the related parame- 
ter which they denote by I c (the overall critical current 
of their nonuniform thickness junction) is found to have 
a significant dip at the <-> tt transition temperature. 
Remarkably, these transitions were observed in Ref. 16 
even though their samples did not have layers of uniform 
thickness. These two experiments and others show that 
this is an observable and robust phenomenon. 

A similar analysis for the 7 layer case showed the same 
correspondence between AJF(0) and T c . However, the 
larger number of energy crossings at T = lead to a 
closer spacing of energy crossings in Dp, causing sev- 
eral local minima in AJ r (0) to appear as a single broad 
minimum. The result was that multiple dips in T c often 
merged. In larger systems, the existence of a broad local 
minimum in T c may correspond to multiple crossings at 
low temperature. 



F. Comparison with experiment 

Many experimental groups 20-24 have found oscillations 
in T c as a function of F layer thickness. Our calculation 
also finds these oscillations (see Fig. 10). The agreement 
with experiment is furthermore quantitative. We show in 
Fig. 11a direct comparison of our results to the exper- 
imental data for a Nb/Co system. 20 In the experiment, 
the spontaneous magnetization of the Co layer was found 
to depend on its thickness. This means, in our language, 
that / must be taken as a function of Dp for the pur- 
poses of this comparison. To do this, we fitted the spon- 
taneous magnetization reported in Ref. 20 and extracted, 
at each thickness, the value of I from the formula below 
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FIG. 12: Density of states at T x for an SFS trilayer. The 
quantity plotted is the local DOS as defined in Eq. (2.6), 
averaged over an S layer, and normalized to the normal state 
bulk result in S. This is the case shown in the bottom left 
panel of Fig. 1 (T x /T° = 0.16). 



FIG. 13: The dimensionless local magnetization M(Z) for 
an SFS trilayer at T x . Only the central part of the sample is 
plotted. Results are given for the two nearly coexisting states. 
The case plotted is the same as in Fig. 12. The vertical dotted 
lines delimit the F region. 



Eq. (2.8). We took a constant A = 0.60, which is appro- 
priate to the materials mentioned. The experimental and 
theoretical values are in excellent quantitative agreement 
on the vertical scale, and the damped oscillations are are 
well aligned in the thickness. 

Comparing figures 11 and 10, we conclude that the dips 
in Fig. 11 must correspond to changes in the stable state 
at zero temperature. As these changes are, as we have 
seen, associated with the first order phase transitions, 
these dips in T c may also be associated with first order 
phase transitions, in good agreement with what was re- 
ported in Rcf. 24. This implies that studies of T c may 
be a useful tool for experimental discovery of first order 
phase transitions and that samples which show dips in 
T c are the ones that should be cooled down and studied 
to locate such phase transitions. 



G. DOS and M(z) 

Advanced tunneling spectroscopy techniques are a use- 
ful experimental tool to measure the local DOS, thus 
probing the single-particle spectrum. It has been found 6 
previously that the local DOS results for and ir states 
are different, including a modified subgap structure. In 
such cases, tunneling spectroscopy could be used to dis- 
tinguish the states. We now investigate whether the den- 
sity of states is also a suitable technique in locating phase 
transitions. 

In Fig. 12 we show the DOS, defined as the normal- 
ized local DOS (from Eq. 2.6) averaged over one of the 
S layers, for a typical 3 layer system at the temperature 
where the first order transition occurs. The case shown 
is for the same parameters as in the bottom left panel of 
Fig. 1, with T x /T® — 0.16. The energy is normalized to 
the bulk S gap at zero temperature, Ao, while the DOS 



is normalized to its value in a bulk sample of S material 
in its normal state. For both and n states maxima ex- 
ist near the bulk gap edge, qualitatively reminiscent of 
the divergence found in a bulk superconductor. The local 
DOS never quite goes to zero in either state, demonstrat- 
ing gapless superconductivity induced by the numerous 
Andreev bound states in the gap. The number of states 
in the gap is clearly larger for the state and the peak is 
markedly lower. Although the DOS for both states have 
some general similarities, the differences that do exist are 
well within the resolution of current tunneling probes, 31 
making the DOS a potentially useful experimental tech- 
nique in locating the phase transitions or identifying the 
stable state in the neighborhood of a transition. 

The last quantity we shall briefly describe is the lo- 
cal dimensionless magnetization, M(z), as defined in 
Eq. (2.8). Previous studies 6 indicate that there is lit- 
tle difference between M(z) for the and ir states at low 
temperature. In that case it was found that M(z) was 
dominated by the exchange parameter / and was rather 
insensitive to the phase of the superconducting state. 
There was also little magnetization induced in the S re- 
gion, as M(z) decayed over the Fermi length scale. 6 To il- 
lustrate the effect that temperature has on this trend, we 
show M(Z) versus the dimensionless length Z at T = T x 
in Fig. 13. In the figure, the F region is delimited by 
vertical dotted lines and only a small portion of the S 
regions is shown. Consistent with Ref. 6, there is a quick 
decay and oscillation of M(Z) in the S region. There is 
a rise in the value of M(Z) to about 0.33 in the center of 
the F region, which is consistent with the bulk formula 
below Eq. (2.8) for / = 0.2. Indeed, as Uf increases 
M(Z) flattens to a value that is in good agreement with 
that estimate. This is not contrary to the experimental 
results in Ref. 20 for which we modeled the change in the 
saturation magnetization with Dp by allowing for a Dp 
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dependent /, since in that case the magnetic properties 
(such the saturation magnetization) of the F layer were 
experimentally found to depend on Dp. Thus, the lo- 
cal magnetization, while interesting for other reasons, is 
not a good tool for determining the thcrmodynamically 
stable state or locating phase transitions. 

IV. CONCLUSION 

We have rigorously studied the thermodynamics of 
clean SFS trilayer junctions through self-consistent solu- 
tions to the BdG equations, in the clean limit. Building 
on previous work 35 where <-» n transitions in this sys- 
tem were found to be possible, we have computed here the 
three dimensional phase diagram of a clean SFS junction 
over an extended and physically relevant region of the 
space spanned by the parameters T, dp, and A. We have 
found that the transition to the normal state is always of 
second order, while first order tt — > transitions occur, as 
temperature increases, over a range of A and Dp. Such 
transitions have been found experimentally. For systems 
consisting of three such junctions, we have found here 
that a variety of first order transitions, involving <-> 7r 
switching of one or more junctions, occur. The phase 
transitions were shown to be driven by a delicate bal- 
ance between the condensation energy and the entropy. 
The absolute value of the pair amplitude is discontinu- 



ous at the first order transition. Key elements of our 
approach are an efficient method to accurately compute 
free energies and a linearization scheme that calculates 
T c . We have shown that dips in T c overlap with regions 
in parameter space where phase transitions exist, which 
suggest that T c studies should be useful for experimen- 
tally locating first order phase transitions. We have also 
calculated the variation of T c with dp and found good 
quantitative agreement with an experimental study 20 of 
a Nb/Co system. We have demonstrated that the phase 
transitions will have measurable latent heats, even for 
relatively small samples, over a broad range of magnet 
thicknesses. Another experimentally relevant quantity, 
the DOS, was calculated and deemed a potentially useful 
tool in locating phase transitions. The local magnetiza- 
tion, however, shows little difference between two states 
at the first order transition. The method and results 
demonstrated here are expected to be applicable to even 
larger structures. 
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